GLHT Method 1: Symbolic Description

To compare factor levels, use the mcp() function.

islandTest <- glht(reg4, linfct = mcp(island = "Dream - Torgersen = 0"))
summary(islandTest)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lm(formula = body_mass_g ~ bill_length_mm + flipper_length_mm + 
    island, data = penguins)

Linear Hypotheses:
                       Estimate Std. Error t value Pr(>|t|)  
Dream - Torgersen == 0  -142.76      69.72  -2.047   0.0414 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

GLHT Method 2: Contrast Matrix

We could also specify a contrast vector or matrix:

coef(reg4)

      (Intercept)    bill_length_mm flipper_length_mm       islandDream 
      -4355.02693          16.88279          39.65617        -333.71627 
  islandTorgersen 
       -190.95953 

K <- rbind("Dream - Torgersen" = c(0, 0, 0, 1, -1))
K

                  [,1] [,2] [,3] [,4] [,5]
Dream - Torgersen    0    0    0    1   -1

islandTest2 <- glht(reg4, linfct = K)
summary(islandTest2)

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = body_mass_g ~ bill_length_mm + flipper_length_mm + 
    island, data = penguins)

Linear Hypotheses:
                       Estimate Std. Error t value Pr(>|t|)  
Dream - Torgersen == 0  -142.76      69.72  -2.047   0.0414 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Linearity

With multiple covariates, the linearity assumption becomes even more critical. Not only do we have to correctly model the form of the relationship between Y and X_1, but also X_2, X_3, etc. We can use a plot of residuals vs. fitted values to check this assumption. We could also make partial regression plots to check.

library(car)
crPlots(reg2)

Outliers: Leverage

Leverage is a measure of how much the data from the ith observation contribute to the fitted value \hat{Y}_i. High leverage means that X_i strongly determines \hat{Y}_i. Remember \hat{Y} = HY = {X}^\top \left({X}^\top {X}\right)^{-1} {X} {Y}.

The leverage for the ith observation is the (i,i)th element (the ith diagonal) of the hat matrix H, h_{ii} = x_{i}^\top (X^\top X)^{-1} x_{i}, and reflects outlyingness in X space.

Outliers: Studentized Residuals

The standardized residuals are scaled by an estimate of the errors’ variance. We could also scale by an estimate of the residuals’ variance. The

Remember,

\operatorname{Var}(\hat{\epsilon}) = \sigma^2(I_n - H),

where H = X^\top (X^\top X)^{-1} X is the hat matrix and I_n is the n\times n identity matrix. The diagonal elements of the hat matrix, denoted h_{ii} are called leverage. The variance of a single residual is \operatorname{Var}(\hat{\epsilon}_i) = \sigma^2 (1 - h_{ii}). The studentized residuals are

\hat{t}_{i} = \frac{\hat{\epsilon}_i}{\hat{\sigma}^2 \sqrt{1-h_{ii}}}.

The studentized residuals have constant variance 1 regardless of the location of x_i when the model is correctly specified. Studentized residuals are useful for examining outliers because they’re leverage-adjusted: they account for the fact that the variance of the residuals changes with X_i.

Studentized residuals are approximately distributed with a standard normal distribution: “large” values (per a standard normal distribution) indicate outlyingness.

Outliers: Jackknife Residuals

The jackknife residuals are similar to studentized residuals, but refit to remove the ith individual:

\hat{r}_{(-i)} = \frac{\hat{\epsilon}_i}{\hat{\sigma}_{(-i)} \sqrt{1 - h_{ii}}},

where \hat{\sigma}_{(-i)}^2 is an estimator of \sigma^2 with the ith individual removed (i.e., that doesn’t use data from individual i).

Jackknife residuals completely remove the influence of the ith individual in determining whether it’s an outlier, because \hat{\epsilon}_i can mask its own outlyingness:

• For points with high leverage, |\hat{\epsilon}_i| may be small (remember the variance of the residual is \sigma^2 (I - H))

• If |\hat{\epsilon}_i| is large, \hat{\sigma}^2 might over-estimate \sigma^2, so things like the studentized residuals might be too small.

Influence

Influence measures the “pull” a single point has on a fitted regression line. Highly influential points draw the fitted line closer to them, possibly biasing results. They typically have both

• high leverage

• departure from trend in other data points

but one or the other does not necessarily imply influence. An influential point is one with high leverage and outlyingness.

One measure of influence is Cook’s Distance:

D_i = \frac{\left(\hat{\beta}_{(-i)}-\hat{\beta}\right)^\top \left(X^\top X\right)\left(\hat{\beta}_{(-i)} - \hat{\beta}\right)}{(p+1)\hat{\sigma}^2}.

Cook’s distance is a scaled distance between the fitted regression coefficients after removing observation i and the fitted regression coefficients that use observation i. If D_i > 1, generally we’ll say the observation is influential.

We can plot Cook’s distance for all points in the model using the which argument to plot.lm():

plot(reg2, which = 4)